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Abstract 

The classical level set method, which represents the boundary of the unknown 
geometry as the zero-level set of a function, has been shown to be very effective in 
solving shape optimization problems. The present work addresses the issue of using a 
level set representation when there are simple geometrical and topological constraints. 
We propose a logarithmic barrier penalty which acts to enforce the constraints, leading 
to an approximate solution to shape design problems. 

1 Introduction 

The level set method El E] is a very powerful approach for problems involving geometry 
and geometric evolution. It has also been applied to solving shape optimization problems 
[H HE IH] j and it is at this type of problems that this work is aimed. 

By a shape we mean a bounded region D in M. n with C 1 boundary. One associates with 
D a function (f> : M. n — > R with the property that D is the level set of </>, 

D = {x : <f>(x) > 0}. 

One then manipulates D implicitly, through its level set function <fi. It is typical in shape 
optimization problems to start with an initial shape, which is then improved in an iterative 
process. Thus, one would start with a level set function 0(x) which is updated at each 
iteration. 

The advantage of the level set method is that it is much easier to work with a globally 
defined function than to keep track of the boundary of a domain. The latter, which can be 
achieved by using marker points and spline interpolation, can become especially complicated 
if D has either several connected components, or is otherwise connected but has several 
holes. During the optimization process, the components or holes may merge or split, or 
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even entirely disappear. The level set method, on the other hand, takes care of this kind of 
changes with great ease. 

Given the shape D there exist of course many functions whose level set is D. The most 
convenient <p to work with is the signed distance to the boundary dD of D, thus 



for x in a neighborhood of dD. Any level set function </> can be reinitialized as the signed 
distance to the set {x : <p(x) = 0}, so from here on we will assume that (f> always satisfies 
(|T|), by reinitializing it if necessary. 

It is very easy to describe deformations of D in terms of its level set function (p. For 
example, if h : M. n — ► R is a function with sup |/i(x)| small enough, then the level set of + h 
is obtained from the level set D of <fi by shifting every point x € <9.D by approximately the 
amount h(x) in the direction of the external normal to dD at x (which is — V0(x)). 

While the level set method has its strong points - one being that it gives a representation 
that is topology-independent - it is not obvious how to extend it to problems where there 
are constraints. Simple volume (area in 2-D) constraints are relatively easy to incorporate 
jSj- Other constraints, such as a bound on the size of a connected component of D, or the 
requirement that D has a fixed number of connected components, are not as easy to handle. 
It is towards this class of problems that this work is directed. 

Our approach starts with the concept of subdomain neighborhood. The neighborhood 
of one subdomain will detect the nearness of other subdomains, and will thus allow us to 
take action to prevent geometry or topology changes. This strategy can be formulated as a 
penalty functional, which we describe in the next section. We illustrate this method by two 
numerical examples in Section 3. 

We wish to mention the paper [2] which also suggests a way of adapting the level set 
method to preserve topology. The authors of this paper do it in the context of image 
segmentation. The key difference between our work and |2] is that their method is pixel- 
based. The algorithm in [2] is able to detect that a shape is about to change topology only 
when certain dimensions of the shape are of size comparable to the grid size. In the context 
of image processing this makes a lot of sense, as then it is convenient to define a body to 
be connected as long as it is made up of one or more pieces joined together by at least one 
pixel. 

We developed our topology preserving level set method having in view problems of shape 
design. There, one specifies in advance certain conditions about how small, thin, or close 
certain features of the shape can get, and then one uses a grid as fine as needed to resolve 
the details of the optimal shape. Thus, our method will be different from j2] by the fact that 
our method is grid size independent. 




x e D, 

x^D. 



(1) 



Then <p will have the additional property 



V<f>{x) ■ V<f>(x) 
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(2) 
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2 Topology-preserving level set method 



A typical shape optimization problem is as follows. We are given a cost function F which 
depends on geometry of the unknown shape. The problem is to find a shape such that the 
cost function is minimized (at least locally). 
Let us represent the shape D as 

D = {x : 4>(x) > 0}. 

The optimization problem we wish to solve is 

min F ((/)), 
4> 

subject to geometrical and topological constraints on D. The latter constraints are: 

• Shape topology. The domain we design for must have, for example, a fixed number 
of connected components or holes. 

• Component size. A lower bound on the size of each component or hole is prescribed. 

• Distance between components. A lower bound on the distance between compo- 
nents or holes is prescribed. In the case of holes, we also prescribe a lower bound on 
the distance from each hole to the external boundary of the domain. 

These constraints arise naturally in optimal design problems as we will illustrate in two 
numerical examples. 

It turns out that all these constraints can be handled in a single penalty formulation. We 
will restrict our attention to 2-D problems, even though the same ideas will work in higher 
dimensions. 

Assume for simplicity that D is a bounded and connected set in IR 2 with a set of holes 
inside of it, which are connected components of M. 2 \D. If d > and I > are real numbers, 
denote 

I d = {x + dV<p(x) : x E dD}, 

and 

Ei = {x - lV(f){x) : x E dD}. 

It follows from (J2J) that for d and I small enough, I d and E\ are made up of points at distance 
d and I respectively from dD. In fact, for d = I, the union of these two sets is exactly the 
set of all points at distance d from dD. Note that if any two components of M. 2 \D (we 
consider the unbounded component too) are at distance more than d from each other, then 
Id is entirely inside of D, and thus 4>{x) > on Id- Also, if the gaps in D are not too "small" 
or too "thin", then Ei is a subset of M?\D, and so <p(x) < on E\. 

Then, we claim, and using a little bit of geometric intuition it is easy to see that it is so, 
that for d and / small numbers, the conditions 

<j)(x + dV(j)(x)) > and 4>{x - ZV^(x)) < for x E dD 
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Figure 1: The set Id (dashed curves) and E\ (dotted curves). 



are a reasonable way of guaranteeing that the holes in D will not merge, split, or become too 
small. In fact, since the outer boundary of D is defined by the same level set function, the 
above also ensures that the holes will never get too close to the boundary. These properties 
also guarantee that if we start an iterative process with the desired topology, the iterations 
cannot change the topology of D as it is updated. Thus, these two conditions on the level 
set function cj) achieve the constraints of the problem. 

To incorporate these conditions into the optimization problem we use the logarithmic bar- 
rier method, see Instead of trying to minimize F((p), consider the problem of minimizing 
F e (</>) = F(<f>) + eH(<j>) for e < 1, where 

H((f>) = - y log [<p{x + dV(j){x))] ds- J log [ - (p(x - lV(f){x))] ds. 

dD dD 

To obtain cj) minimizing F e (<f)) we will use the steepest descent method. It amounts to 
finding the derivative of F £ (<j)), and at each iteration taking a step in the direction in which 
the derivative decreases fastest. 

In order to calculate the derivative of F £ ((p) we need the derivatives of F((p) and H((p). 
Let h : M 2 — > M. be a test function. For t a real number, \t\ 1, F((p + th) will depend on 
the values of h only close to the boundary of D, as F is a function of the level set of + th, 
and the way this level set depends on h was discussed above. We deduce that 

dF(<p + th) 



D^F(cP) ■ h 



dt 



t=o 

will only be a function of the restriction of h to dD. In many important applications, see 
[Hj, it has the form 

D^F{4>) ■ h = J U{x)h{x) ds, (3) 

dD 
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dD 



for some function U which of course depends on cj) and which can be calculated numerically. 

The derivative of H(<p) can be calculated explicitly. Consider a parameterization x(s) 
of dD, with x'(s) having unit norm for all s. H(<ft + th) will be a sum of two integrals 
over the set {x : (0 + th)(x) = 0}, which, if (J2J) holds, is approximately parameterized by 
x — th(x)V 4>{x) , with x = x(s). One can then find that the derivative of the first integral in 
H((f> + th) at t = is 

/ [V0(a;)ft(x) + rf/i(x)V 2 0(x)V0(x) - dV/i(z)] • V<p(x + dV0(x)) - h{x + dV(f>(x)) 
\ <j)(x + dV(p{x)) 

+ log [<fi(x + dV<p(x))]x ■ [{Vh(x) ■ x')V(p{x) + (V 2 (f)(x)x') h(x)] J ds. (4) 

A similar equality holds for the second term in H((p). 

Beside the obvious complexity of this expression, note that unlike the case of F((f>), this 
derivative will no longer depend on the values of the test function h only on dD. We will 
make several approximations. Recall that the purpose of H((p) is to make sure at every step 
in the optimization process the domain D has the topology preserved. H((f>) will grow large 
only when D is close to violating the restrictions imposed on it. As far as the first integral 
in H((f>) is concerned, this happens when <p( x + dS7(f)(x)) becomes close to zero. Then, the 
term on the first line of (j3J) is much larger than the second. We will ignore the term on the 
second line. Also, on the first line, we have V 2 0(x)V0(x) = 0, which follows from (J2J). In 
addition, we will ignore the quantities — dVh(x) ■ V(j)(x + dSJ(j){x)j and —h(x + dS7(j)(x)). 
We obtain the more manageable expression 



j Ui(x)h(x)ds, 



dD 

with 

, x V0O) • V^(x + rfV0(x)) 

The derivative of the second integral in H can be calculated, and then approximated, in 
the same way. Make the notation 

u i i v^) ■ vjK; - 'W(^)) c „ n ... 
' Mx) = — 0(»-tv^» — ' x e aD - (6) 



We obtain 



D^H(<P) -h = J [Ui(x) + U 2 (x)]h(x) ds. 



dD 

This gives us the following approximate equality 



D F £ (0) • h = j [U(x) + eE7i(x) + eU 2 {x)] h(x) ds. 



dD 
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If this were an exact equality, the steepest descent direction for F E (4>) at would be 

u(x) = -\U{x) + eU^x) + eU 2 (x)] , (7) 

where x G 3D. This quantity can be extended continuously to a neighborhood of 3D in 
the following manner: for x G M 2 close to dD let x G dD be the unique point such that 
dist(x,<9D) = dist(x,af), and set 

u(x) = 4>{x) + u(x). (8) 

Then the next iteration for would be + au, where a > is the length of the step to be 
taken in the direction u. 

But the obtained u is an approximation. It will then clearly not be the steepest descent 
direction for P e (0). One could question if it would be a descent direction at all, that is, 
whether P e (0) would decrease if is replaced by + au. After a numerical study we can say 
that the answer is no; P e (0) could even increase in the process. Nevertheless, we will argue 
below that this iterative process does its job at maintaining the topology constraints. And 
as far as the problem of minimizing P(0), it is clear that the iterative process we suggest 
will give us a sufficiently good approximation to the point of minimization 0, provided that 
e is small enough. 

We will show that, if the level set function is such that two components of {x : 4>{x) < 0} 
are at distance slightly more than d from one another, then u will act as a repelling force, 
and in consequence, the components of {x : (4> + au)(x) < 0} will be further apart. 




Figure 2: A blow-up of Fig. [T] 

Indeed, consider such a situation in Fig. |21 Let x G dCi be a point, which in this 
figure we will denote by Pi, such that dist(Pi, <9C 2 ) is slightly larger than d. Let P 2 be the 
point x + dSJ(j)(x). Then P 2 will be very close to dC2- We will have 4>(x + <iV0(x)) > 
but very small. It is easy to show, and geometrically clear, that the gradient of at P 2 , 
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V0(x + dS7(f)(x)), which in the figure is represented by the vector P2P3, will point almost 

in the opposite direction of -P1-P2 which is dV(f)(x). We find that Ui(x) will be negative and 
large in magnitude. 

Moreover, when the distance between dC\ and <9C*2 is close enough to d, eUi(x) will be 
larger in magnitude than U(x) + eU2(x). In consequence, u(x) defined by © will be positive. 
Therefore, we have 4>(x) = 0, but (0 + au)(x) > 0. The same reasoning applies for points 
x G <9C*2 close to dC\. This shows that the connected components of (<fi + au)(x) < will 
be further apart. 

It can be argued in the same manner that should a component of {x : <p(x) < 0} get too 
"thin" or too "small", then L^(x) will serve as a counterweight, forcing it to get "fatter". 

Let us note that in order for the above to work, each step size should not be too big. If 
the boundary of D moves by more than ^ at some step, then two components of {x : <p(x) < 
0} which were at distance slightly more than d can end up merging without the penalty 
functional noticing that. Or, if the boundary moves by more than |, a connected component 
slightly thinner or larger than I might end up splitting or disappearing. Therefore, at each 
step one needs to make sure that 

a max \u(x)\ < Kmin(d, I) (9) 

xGdD 

for some K > 0, as the quantity on the left determines by how much the boundary of D gets 
shifted at the given step. Theoretically K can be allowed to be as large as |, but since we 
use a finite grid size we have to be more conservative. A value of K = | works in practice. 

But enforcing © is not enough to guarantee our geometrical and topological constraints. 
The penalty functional H((f>) is supposed to take care of this, but it is clear that the smaller 
e is, the weaker the influence of H(<f)) in F £ ((p) will be, and the closer to violating the 
constraints <p will get, before this penalty functional kicks in. Thus, at each iteration one 
needs to first take a step size a satisfying ©, and still check after updating <fi to <p + au 
whether H(<fr) is defined. If not, one needs to decrease the step size a, for example by halving 
it, until H((p) is defined. If no amount of decreasing a helps, one needs to either increase e 
or decrease the grid size, and restart the algorithm. 

A pseudo-code for the algorithm is as follows. 



initial guess for <fi 
do while not optimal 

• compute the descent direction u (use (jSJ) , © , © > © > an d ©) 

• choose a step size a satisfying © for which H(<p + au) is defined 

• update (f) to + au 

• reinitialize to satisfy © 



We note that if at some point the contour {x : 4>(x) = 0} develops sharp angles, then the 
functional H(<p) might not be defined (this can be seen from Fig. ^). To prevent this from 
happening, one can smooth <fi a bit at each iteration. For cf) discretized on a square grid we 
used the procedure 

, <f>i,j + <Pi-i,j + <Pi+i,j + <kj-i + <Pi,j+i 
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Figure 3: The initial and optimized shape for example 1. 



Also, for fine grids it becomes expensive to reinitialize <j) according to (JTj). To make this 
computation faster we reinitialized 4> only in a neighborhood of the set {x : 4>{x) = 0}. For 
more performance one could use the fast re-distancing algorithms suggested in [7[ ITUl ITT]. 

Lastly, sometimes one might wish to introduce additional constraints of the form G((p) = 
const, in the optimization problem. An example of such a constraint is the requirement that 
the area of the set {x : (p(x) > 0} be kept fixed, which we will use in the two numerical 
examples below. Then one needs to modify the descent direction u as described in 



3 Numerical examples 

In the first example, we consider the problem of finding a domain that has the smallest 
perimeter, subject to the constraint that the area of the domain being fixed. Thus, the 
functional to minimize is 




with the constraint 

G{4>) — J 1 dx dy = const. 

{0>O} 

The starting shape is a region with seven subdomains, each one an ellipse with aspect ratio 
1.3, as shown in Fig. El on the left. The center ellipse has a slightly bigger (20%) size than 
the rest. The distance between the centers of the ellipses is 4, and the smallest semi-axis of 
the surrounding ellipses is 1. 

If we do not constrain the topology or geometry, the optimal solution would be a circle 
whose area is equal to the area of the original seven subdomains. If we do enforce these 
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Figure 4: The initial and optimized shape for example 2. 



constraints, minimizing instead the functional 

we obtain the picture in Fig. El on the right. 

For this calculation we set d = I = 0.8, e = 0.2 and consider a square grid of size 
h = 0.05 (each square is further split into two triangles, to make it easier to keep track of 
the set {x : (f>(x) = 0}). 

We find that the "satellite" components of the central domain do not disappear, but 
became of size slightly larger than I. 

We note that that the resulting large domain in the center is not perfectly circular. This 
because the steepest descent direction for F(cf>) will be Acf) = <p xx + <p yy . We need to calculate 
this quantity numerically, and after reinitializing <fi according to ((TJ) it is not smooth enough 
for Acj) to be calculated accurately. Smoothing (p as noted in the previous section helped a 
bit, this is how this picture was obtained. We found that if we perform additional smoothing 
then the result in Fig. El will look more circular. This artifact does not show up in the next 
example, as then one does not need to calculate second-order derivatives of <fi. 

In the second example we examine the problem of minimizing the functional 



F{4>) 



[x 2 + y 2 ) 



dx dy. 



{4»o} 



We again enforce the area constraint G{4>) = const., and we use the same values for d, I and 
h. We set e = 0.4. (The value of e which is relatively small, and in the same time be not 
small enough that the algorithm fails converge for a given grid size is determined by trial 
and error, and it depends on the problem.) In absence of topological constraints, these seven 
ellipses would merge to form a large circle. The topological constraints prevent them from 
doing so, as we see from Fig. 0] 



9 



4 Discussion 



In this paper we introduced a penalty functional which makes it possible to use the level 
set method in problems with topology and geometry constraints. Our method allows for 
topological constraints independent of the grid size (that is, for given d and I, the grid size 
h can be chosen as small as desired), which is a key difference with the method suggested in 

m 
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